
import commands
import sys
from os.path import isfile
from scipy import array,isnan,ones,arange,unique,shape,nan,concatenate,diff,sign,transpose,argmax,argmin,argsort,sum,polyfit
from matplotlib.mlab import find
from pylab import figure, axes,xlabel,ylabel,yscale,plot,ylim,xlim,twinx

import IV
from loadDATfile import loadDATfile

class IVTUdelft(IV.IV):
  """
  Do analysis of IV data from TU delft fileformats. 
  
  """
  def loadMCBJ(self,basename,force=False):
    """
     Load IV files from dir rawdata/"basename"*.dat
    """
    self.basename=basename[len('rawdata/')+basename.find('rawdata/'):]
    if isfile("%strace.dat" % (basename)):
      self.trace = loadDATfile("%strace.dat" % (basename))
      self.SavedTracestmp=self.trace[1][self.trace[1][:,3]<>-1]
    else:
      print "Trace file is missing: %s" % ("%strace.dat" % (basename))
      self.trace = None
      self.SavedTracestmp  = [None]
    if isfile(basename+'.npz') and (not force):
      self.loadz(basename+'.npz')
      print "Import from %s. To import from dat files use force=True" % (basename+'.npz')
      return True
    filer=commands.getstatusoutput("ls -t -r "+str(basename)+"*.dat")[1].split("\n")
    Nfiles = len(filer)
    print "Detected %i IV's." % ((Nfiles-1)/2.)
    data=[]
    self.motor=[]
    self.meta = None
    SavedTraces=[]
    for i in range(0,Nfiles):
      fil = str(basename) +str(i)+ ".dat"
      if not isfile(fil):
        continue
      if not isnan(self.SavedTracestmp).any(): 
        SavedTraces.append(self.SavedTracestmp[self.SavedTracestmp[:,3]==i,:][0] if  sum(self.SavedTracestmp[:,3]==i)==1 else [nan,nan,nan,i,nan,nan])
      if ((i%(1+len(filer)/100)) == 0):
        sys.stdout.write('+')
        sys.stdout.flush()
      meta,dat = loadDATfile(fil)
      if len(dat.shape)==2:
        data.append(dat)
        self.motor.append(float(meta['Axis position']))
        if self.meta is None:
          self.meta = meta
    self.motor = array(self.motor)
    sd = array(map(shape,data))
    sdatdict = dict(sd)
    if len(sdatdict)>1 or sum(sd.std(0)>0):
      shapescount = array([(sum((sd[:,0]==i)*(sd[:,1]==sdatdict[i])),i,sdatdict[i]) for i in sdatdict])
      maxshape = shapescount[argmax(shapescount[:,0]),1:]
      print "dataset contains multiple shapes of datafiles:"
      print sdatdict
      print "Only the shape with most files is sellected"
      print maxshape
      self.data  = filter(lambda dat: (shape(dat)[0]==maxshape[0]) & (shape(dat)[1]==maxshape[1]), self.data)
      self.motor = self.motor [(sd[:,0]==maxshape[0])*(sd[:,1]==maxshape[1])]
    data = array(data)
    self.SavedTraces = ones((len(data),6))
    self.SavedTraces[:,0] = arange((len(data)))
    self.SavedTraces[:,3] = arange((len(data)))
    SavedTraces=array(SavedTraces)
    doS = unique(concatenate((find(diff(sign(diff(self.motor)))<0)[::-1],SavedTraces[:,3][diff(SavedTraces[:,4])>0])))
    for i in doS:
      self.SavedTraces[1+i:,4] +=1
    self.SavedTraces = array(self.SavedTraces)
    self.V=data[0,:,0]
    for j in range(2):
      self.Imean[j][0] = data[:,:,j*len(data[0,0,:])/2+1]
      self.I[j][0] = transpose(data[:,:,j*len(data[0,0,:])/2+2:(j+1)*len(data[0,0,:])/2],(0,2,1))
    self._doNs_()
    self.savez(basename+'.npz')
  def PlotIM(self,traces=None,plotlowbias=True,plotmotor=True,**kw):
    figure(figsize=[ 12.90,   7  ])
    ax1 = axes([0.1,0.33,0.8,0.6]) if plotlowbias else axes([0.1,0.15,0.8,0.8])
    IV.IV.PlotIM(self,traces=traces,**kw)
    Xlim=xlim()
    if plotlowbias:
      ax1.xaxis.set_ticks_position ('top')
      ax2 = axes([0.1,0.1,0.8,0.2],sharex=ax1);xlabel('Trace nb.')
      if not vars(self).has_key('Gfit'):
        self.zerobiasconductance()
      plot(array(traces),self.Gfit[traces],'-r') if traces is not None else plot(self.Gfit,'-r')
      yscale('log');ylabel(r'low bias $\frac{G}{G_0}$')
      YL=ylim()
      for i in find(diff(self.SavedTraces[traces if traces is not None else range(len(self.SavedTraces)),4])>0)+.5:
        plot([i,i],YL,'k--',lw=2)
      ylim(YL)
    if plotmotor:
      ax3 = twinx()
      plot(array(traces)+.5*0,self.motor[traces],'--b')  if traces is not None else plot(self.motor,'--b') 
      t = ax3.yaxis.get_ticklocs()
      ax3.yaxis.set_ticks(t[::len(t)/4])
      ylabel(r'Motor position')
      xlim(Xlim)
  def _doNs_(self):
    N = arange(len(self.V))
    m1 = N[argmax(self.V):argmin(self.V)+1*sign(argmin(self.V)-argmax(self.V)):sign(argmin(self.V)-argmax(self.V))]
    m2 = filter(lambda n: not any(m1==n),N)
    M = m1,array(m2)
    M = [m1[argsort(self.V[m1])] for m1 in M ]
    if self.V[M[0][0]]<self.V[M[1][0]]:
      M[1] = array([M[0][0]] +list(M[1]))
    if self.V[M[0][-1]]>self.V[M[1][-1]]:
      M[1] = array(list(M[1]) + list([M[0][-1]]))
    M = [m1[diff(self.V[m1])<>0] for m1 in M ]
    self.Ns = M
 